import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint, solve_ivp

def ecosystem_model(t, y, params):
    I, P, A, B, Bee = y
    r_I, K_I, r_P, K_P, r_B, K_B, r_Bee, K_Bee, beta_P, alpha_P, alpha_A, tau_B, eta_B, sigma_B, zeta_B, mu_B = params

    dI_dt = r_I * I * (1 - I / K_I) + beta_P * P - alpha_P * P * I - alpha_A * P * A - eta_B * I - zeta_B * I * B
    dP_dt = r_P * P * (1 - P / K_P) - alpha_P * P * I - alpha_A * P * A + tau_B * B * P
    dA_dt = 0
    dB_dt = r_B * B * (1 - B / K_B) + eta_B * I - mu_B * B
    dBee_dt = r_Bee * Bee * (1 - Bee / K_Bee) + sigma_B * P

    return [dI_dt, dP_dt, dA_dt, dB_dt, dBee_dt]

def objective(params):
    I0 = 50
    P0 = 100
    A0 = 30
    B0 = 5
    Bee0 = 5
    y0 = [I0, P0, A0, B0, Bee0]

    t = np.linspace(0, 50, 500)

    r_I, K_I = 0.1, 1000
    r_P, K_P = 0.2, 500
    r_B, K_B = 0.1, 100
    r_Bee, K_Bee = 0.1, 50
    beta_P, alpha_P, alpha_A = 0.05, 0.02, 0.01
    tau_B, eta_B, sigma_B = params[0], params[1], params[2]
    zeta_B, mu_B = 0.01, 0.1

    model_params = [r_I, K_I, r_P, K_P, r_B, K_B, r_Bee, K_Bee, beta_P, alpha_P, alpha_A, tau_B, eta_B, sigma_B, zeta_B, mu_B]
    result = solve_ivp(ecosystem_model, [0, 50], y0, args=(model_params,), t_eval=t, method='RK45', atol=1e-8, rtol=1e-8)

    P_final = result.y[1, -1]
    return -P_final  # 我们希望最大化 P_final，因此返回其负值进行最小化

# 优化参数 tau_B, eta_B, 和 sigma_B
from scipy.optimize import minimize

initial_guess = [0.1, 0.1, 0.1]
bounds = [(0, 1), (0, 1), (0, 1)]
result = minimize(objective, initial_guess, bounds=bounds)

# 提取优化后的参数
optimal_params = result.x
tau_B, eta_B, sigma_B = optimal_params

# 使用优化后的参数模拟生态系统模型
r_I, K_I = 0.1, 1000
r_P, K_P = 0.2, 500
r_B, K_B = 0.1, 100
r_Bee, K_Bee = 0.1, 50
beta_P, alpha_P, alpha_A = 0.05, 0.02, 0.01
zeta_B, mu_B = 0.01, 0.1
model_params = [r_I, K_I, r_P, K_P, r_B, K_B, r_Bee, K_Bee, beta_P, alpha_P, alpha_A, tau_B, eta_B, sigma_B, zeta_B, mu_B]
y0 = [50, 100, 30, 5, 5]
t = np.linspace(0, 50, 500)
result = solve_ivp(ecosystem_model, [0, 50], y0, args=(model_params,), t_eval=t, method='RK45', atol=1e-8, rtol=1e-8)

# 绘制结果
plt.figure(figsize=(10, 6))
plt.plot(t, result.y[0], label='Insects (I)', color='blue')
plt.plot(t, result.y[1], label='Plants (P)', color='green')
plt.plot(t, result.y[2], label='Herbivores (A)', color='red')
plt.plot(t, result.y[3], label='Bats (B)', color='purple')
plt.plot(t, result.y[4], label='Bees (Bee)', color='orange')
plt.xlabel('Time (Years)')
plt.ylabel('Population')
plt.title('Ecological System Dynamics with Bat and Bee Introduction')
plt.legend()
plt.grid(True)
plt.show()